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Abstract 

An  advantage  of  linear  quadratic  regulator  (LQR)  design  is  that  it  gives  a  robust 
system  by  guaranteeing  stability  margins.  This  property  is  used  to  develop  an  algorithm 
for  placing  robust  poles.  The  algorithm  chooses  the  positive  semidefmite  weighting 
matrix  Q  and  positive  definite  weighting  matrix  R  by  attempting  to  place  closed  loop 
poles  at  a  set  of  desired  poles.  If  the  desired  poles  lie  outside  the  allowable  LQR  region, 
the  algorithm  finds  the  achievable  poles  inside  the  region  that  are  closest  to  the  desired 
poles.  The  solution  requires  using  a  gradient  search  technique  to  minimize  a  weighted 
eigenvalue  difference  cost  function.  The  weighting  of  the  eigenvalue  difference 
establishes  the  relative  imoortance  between  the  poles.  In  a  multi-input  multi-output 
system,  the  placement  of  one  pole  effects  the  allowable  placement  region  of  the  other 
poles.  Thus,  the  heavier  weighted  poles  have  precedence  and  are  forced  closer  to  their 
desired  location.  The  algorithm  is  programmed  to  run  on  the  software  package  MATLAB 
and  the  related  subroutines  are  discussed.  Several  examples  are  included  to  illustrate  the 
use  of  the  algorithm,  some  of  which  can  be  solved  in  closed  form  to  compare  with  the 
algorithm’s  solution.  The  results  show  that  this  technique  is  accurate  for  selecting  robust 
poles  at  or  close  to  the  desired  pole  locations.  ,  ' 
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A  LINEAR  QUADRATIC  REGULATOR  WEIGHT  SELECTION 


ALGORITHM  FOR  ROBUST  POLE  ASSIGNMENT 


I.  Introduction 

Classical  flight  control  system  (FCS)  design  has  evolved  from  the  need  to  make 
aircraft  stable  and  reduce  pilot  workload.  Early  designs  built  stability  into  the  airframe 
making  them  less  maneuverable  and  more  susceptible  to  atmospheric  disturbances. 
Before  World  War  L  rhe  first  FCS’s  were  simple  gyroscopic  autopilots  which  provided 
stability  and  a  return  to  the  desired  attitude  after  disturbances  f  1 : 1  ].  With  the  increasing 
complexity  of  modem  jets,  modem  FCS's  have  become  an  integral  part  of  the  aircraft. 
Modem  FCS's  not  only  provide  stability  and  desired  responses  to  specific  inputs,  but  also 
suppress  disturbances,  component  variations,  uncertainties,  and  cross  coupling  effects 
[2:27], 

In  the  preliminary  design  phase,  one  of  the  goals  of  the  control  engineer  is  to 
design  the  FCS  to  allow  the  airplane  to  accomplish  its  mission.  The  mission  requirements 
can  be  broken  down  into  several  basic  functions  which  are  quantified  by  specifications. 
Aircraft  specifications  are  usually  governed  by  documents  such  as  Mil  Standard  1797  and 


Federal  Aviation  Regulations. 


The  dynamics  of  the  airplane  can  be  approximated  by  mathematical  modeling. 
The  resulting  differential  equations  establish  the  characteristic  equation  of  the  system. 
To  evaluate  the  system,  the  characteristic  equation  is  factored  into  roots  and  plotted  in  the 
complex  plane.  These  points  are  called  poles  and  they  describe  the  performance  and 
stability  characteristics  of  the  airplane.  The  math  model  of  the  airplane  is  called  the 
plant,  and  certain  variables  or  states  (such  as  roll  rate  or  bank  angle)  of  the  plant  output 
can  be  fed  back  into  the  plant  input  as  shown  in  figure  1 .  Using  feedback,  the  poles  of 
the  system  can  be  shifted.  In  this  fashion,  the  design  engineer  has  the  ability  to  shape  the 
closed  loop  performance  characteristics  of  the  system  to  meet  the  design  goals.  A  variety 
of  FCS  design  techniques  exist;  this  report  is  concerned  with  pole  placement  and  the 
linear  quadratic  regulator  (LQR)  method. 

Background 

Pole  placement  is  a  powerful  tool  for  the  control  system  designer.  In  fact,  a 
controllable  system  can  be  forced  to  satisfy  any  desired  set  of  poles  with  the  appropriate 
linear  feedback  [3].  Thus,  the  designer  can  shape  the  system  characteristics  to  meet  the 
design  goals  by  obtaining  a  feedback  matrix  such  that  the  closed  loop  poles  approach  the 
desired  poles.  In  the  single-input  case,  the  feedback  matrix  is  unique  for  a  given  set  of 
poles.  If  the  poles  are  placed  exactly  at  the  desired  locations,  the  unique  gains  give  fixed 
robustness  characteristics  that  the  designer  must  live  with.  In  the  multi-input  case 
however,  there  is  not  one  unique  feedback  gain  matrix  that  will  yield  a  given  set  of  poles. 
In  other  words,  different  feedback  matrices  will  give  the  same  poles  (or  eigenvalues),  but 
each  one  of  these  feedback  matrices  will  have  a  different  mode  shape  (or  eigenvector)  for 
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Figure  1  State  Feedback  Diagram 


the  system.  For  a  multi-input  system,  Moore  has  shown  that  after  the  poles  are  placed 
an  eigenvector  set  can  be  chosen  from  a  class  of  allowable  eigenvector  sets  [4],  This 
offers  the  freedom  to  first  shift  the  poles  and  then  choose  an  allowable  mode  shape  from 
those  available  A  pole  shifting  technique  for  multi-variable  feedback  given  by  Retallack 
and  MacFarlane  shows  how  to  select  the  feedback  gain  matrix  to  get  any  desired  pole  [5]. 

Even  though  their  pole  placement  method  can  give  the  desired  performance 
characteristics,  it  does  not  guarantee  a  robust  system.  That  is,  one  that  is  insensitive  to 
system  parameter  variations  and  external  disturbances.  Another  way  to  think  of 
robustness  in  a  system  is  that  it  remains  stable  even  in  the  case  of  model  approximation 
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errors  or  flight  condition  disturbances.  A  typical  measure  of  robustness  in  the  frequency 
domain  is  the  gain  and  phase  margin. 

There  are  pole  placement  methods  that  do  achieve  a  robust  system.  However,  the 
poles  are  usually  restricted  to  a  specific  region  of  the  complex  plane.  Simonyi  and  Loh 
developed  a  pole  placement  technique  that  places  the  closed  loop  poles  into  a  region 
where  the  robustness  to  perturbations  is  determined  with  a  condition  number  of  the 
eigenvalue  matrix  [6].  The  linear  quadratic  regulator  (LQR)  is  an  optimal  design 
technique  that  guarantees  a  robust  system.  Anderson  and  Moore  have  shown  that  the 
LQR  design  guarantees  robustness  by  giving  a  gain  margin  of  (-6,°°)  db  and  a  phase 
margin  of  (-60,  60)  degrees  [7:sect  5.4],  The  difficulty  lies  in  choosing  a  weighting 
matrix  for  the  LQR  cost  function  that  gives  the  desired  poles,  and  a  great  deal  of  research 
has  been  accomplished  in  this  area.  Kawasaki  and  Shimemura  [8]  and  Wittenmark  et.  al. 
[9]  have  developed  procedures  to  select  weighting  matrices  to  place  poles  in  a  specified 
region.  Methods  to  place  desired  poles  by  finding  the  state  weighting  matrix  are  given 
by  Graupe  [10],  Broussard  [11],  and  Solheim  [12].  In  papers  by  Wilson  and  Cloutier  [  1 3] 
and  Sobel  and  Shapiro  [14],  a  cost  function  is  minimized  that  has  a  trade-off  between 
desired  pole  placement  and  desired  eigenvector  assignment.  A  method  described  by 
Harvey  and  Stein  [15]  uses  the  state  weighting  matrix  to  place  the  transmission  zeros  in 
the  desired  pole  locations.  As  the  gains  are  increased,  the  closed  loop  poles  migrate 
towards  the  desired  locations.  In  this  report,  an  algorithm  is  developed  that  finds  the 
LQR  state  weighting  matrix  to  place  a  set  of  desired  poles.  The  procedure  varies  the 
weighting  matrices  Q  and  R  to  find  the  smallest  difference  between  an  LQR  achievable 
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pole  and  the  desired  pole.  In  contrast,  Harvey  and  Stein  calculate  Q  to  place  the 
transmission  zeros  and  vary  R  by  defining  R=pl  and  letting  p  tend  towards  zero.  While 
some  of  their  poles  move  to  the  desired  locations,  others  escape  to  infinity.  The 
algorithm  given  here  places  all  the  poles  near  their  desired  locations,  and  none  go  to 
infinity.  Also,  the  algorithm  developed  in  this  report  allows  a  weighting  to  be  placed  on 
selected  poles  to  indicate  their  relative  importance.  A  similar  technique  presented  by 
Graupe  finds  Q  by  minimizing  the  pole  difference,  but  he  does  not  include  any  weighting 
option  or  vary  R.  In  Broussard’s  paper,  an  algorithm  is  presented  to  find  Q  by 
minimizing  the  feedback  gain  difference,  but  the  gains  for  the  desired  poles  must  first  be 
calculated  Broussard  does  introduce  a  weighting  function  in  his  routine,  but  it  is  used 
to  weight  the  eigenvectors  -  not  to  weight  selected  eigenvalues.  The  method  in  this  report 
does  not  calculate  any  eigenvectors.  Instead,  Moore’s  result  stated  above  provides  the 
option  to  select  an  eigenvector  set  after  this  pole  placement  procedure  is  employed.  The 
algorithm  presented  here  is  simple  to  use  but  gives  a  powerful  and  useful  technique  for 
pole  placement. 

Problem  Statement 

Consider  the  state  space  representation  of  a  linear  time  invariant  system 


x  =  Ax  +  Bu 

(1) 

•-c 

II 

n 

X 

(2) 

where  x  is  an  n -dimensional  state  vector,  u  is  an  m-dimensional  input  vector  and  y  is  an 
1-dimensional  output  vector.  A,  B,  and  C  are  constant  matrices  of  the  appropriate  size. 
Assuming  (A,B)  is  controllable,  a  linear  feedback  of  the  state  variables 
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u  =  Kx 


(3) 


can  be  found  that  shifts  the  closed  loop  poles  to  the  desired  locations.  A  problem  arises 
when  a  control  is  required  to  drive  the  plant  output  from  a  nonzero  state  to  the  zero  state 
for  a  plant  subjected  to  unwanted  disturbances.  The  problem  is  to  drive  the  system  back 
to  the  zero  state  as  fast  as  possible  while  trying  to  use  the  smallest  control  inputs. 
Anderson  and  Moore  define  this  as  the  regulator  problem  [7:8].  The  optimal  solution  is 
to  find  the  control  input  u  to  minimize  the  quadratic  cost  function 

J  =  f  (x'Qx  +  u'Ru)  ch  (4) 

where  Q  and  R  are  weighting  matrices  chosen  by  the  designer.  This  is  the  classical  linear 
quadratic  regulator  design  problem.  An  advantage  of  using  LQR  design  is  that  the  system 
will  always  be  stable  and  robust.  The  optimal  solution  is  given  by 

K  =  R  B  P  (5) 

where  P  is  found  from  the  algebraic  Riccati  equation 

PA  +  ATP  +  Q  PBR  BtP  =  0  (6) 

The  LQR  method  guarantees  robustness,  but  only  allows  pole  placement  in  a 
specific  region,  and  the  poles  that  give  the  desired  performance  may  or  may  not  be  in 
these  regions.  Therefore  it  may  not  be  possible  to  use  this  method  to  achieve  both 
robustness  and  exact  pole  placement.  In  that  case  a  choice  would  have  to  be  made.  This 
paper  takes  the  stance  that  robustness  has  higher  priority  over  pole  placement  and  hence 
uses  the  LQR  technique  to  choose  the  poles.  When  the  desired  pole  is  outside  the 
allowable  LQR  region,  a  gradient  search  is  used  to  find  the  achievable  pole  inside  the 
region  that  is  the  closest  to  the  desired  pole.  This  is  done  by  minimizing  the  cost 
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function 


j'  -  E  (7) 

i=l 

where 

•  V,  is  a  positive  definite  weighting  parameter  for  the  poles, 

•  j  is  the  i'h  desired  pole,  and 

^"ach  i  is  the  ilh  achievable  pole  from  LQR 

Thus,  when  J  is  made  small,  the  system  comes  close  to  satisfying  the  desired  pole 
spectrum  and  is  simultaneously  robust. 

Pole  weighting  is  included  in  the  cost  function  to  give  priority  to  poles  that  need 
to  be  a  specific  value.  Some  poles  may  have  limitations  that  prevent  deviations  from  the 
desired  locations.  For  example,  an  actuator  may  have  characteristics  that  determine  the 
location  of  the  pole,  and  if  a  designer  were  to  move  the  pole  he  would  violate  the 
physical  model.  To  circumvent  this,  a  high  weighting  is  placed  on  that  particular 
eigenvalue  difference  in  the  cost  function. 

The  LQR  pole  placement  and  the  pole  difference  minimization  functions  have  been 
compiled  in  an  algorithm  that  finds  the  achievable  poles.  The  algorithm  is  programmed 
into  macros  called  MATLAB  "m"-files  that  run  on  the  software  package  MATLAB.  Five 
test  cases  are  run  with  the  algorithm.  First  and  second  order  systems  that  can  be  solved 
in  closed  form  are  compared  with  the  algorithm  poles;  a  third  order  system  is  run  to  show 
the  effect  of  weighting  an  actuator  pole;  a  sixth  order  example  of  an  F-4  from  Harvey  and 
Stein  is  examined  and  compared  to  their  results;  and  an  example  of  an  aircraft  similar  to 
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an  A-4D  is  examined  to  demonstrate  the  effects  of  robust  pole  placement  on  level  1 
flying  qualities. 

Organization 

The  rest  of  the  report  is  contained  in  chapters  II  through  V  and  appendices  A  and 
B.  They  are  organized  as  follows: 

•  Chapter  II  presents  the  theory  of  the  linear  quadratic  regulator  method  and 

robustness  measures  of  single-input  single-output  (SISO)  and  multi-input 
multi-output  (MIMO)  systems.  A  section  with  definitions  of  observability 
and  controllability  is  included  to  allow  a  complete  presentation  of  LQR 
theory. 

•  Chapter  III  develops  the  robust  pole  placement  algorithm.  A  discussion  on 

how  the  algorithm  is  run  on  the  software  package  MATLAB  is  included. 

•  Chapter  IV  examines  five  test  cases  and  discusses  the  results  of  each.  A 
comparison  of  the  Harvey  and  Stein  method  of  pole  placement  to  the 
algorithm  is  presented  and  the  limitations  of  the  Harvey  and  Stein  method 
are  discussed. 

•  Chapter  V  gives  the  conclusions  and  recommendations  for  further  study. 

•  Appendix  A  supplies  a  listing  of  the  "m"  files  used  to  run  the  algorithm  on 
MATLAB. 

•  Appendix  B  gives  aircraft  data  for  example  five  and  includes  the  formulas 
for  the  longitudinal  equations  of  motion. 
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II.  Theory 


Controllability  and  Observability  [17:  sect  2.2.9] 

Before  the  linear  quadratic  regulator  theory  is  introduced,  some  definitions  that  are 
essential  to  the  theory  are  presented.  A  system  is  completely  controllable  if  any  initial 
state  x,,  can  be  shifted  to  any  final  state  x,  in  a  finite  time.  In  other  words,  a  feedback 
matrix  exists  that  can  shift  a  closed  loop  pole  to  any  desired  location.  Controllability  for 
the  linear  system  given  in  equation  (1)  is  determined  by  forming  the  controllability  matrix 
Mc  =  [B  |  AB  |  A’B  |  ...  |  A-'B]  (8) 

where  n  is  the  number  of  states.  The  system  is  controllable  if  the  rank  of  Mc  is  n.  If  the 
system  is  not  ful1  rank,  n-rank(Mc)  gives  the  number  of  poles  that  can  not  be  moved  by 
feedback. 

A  system  is  said  to  be  completely  observable  if  all  the  states  x  can  be  determined 
from  measuring  the  outputs  y.  If  a  state  is  not  observable,  it  does  not  influence  the 
output  and  cannot  be  identified.  For  the  linear  system  described  by  equations  ( 1 )  and  (2), 
the  observability  matrix 

M„  =  [CT  |  ATCT  |  (AT):CT  |  ...  |  (AT)n  ’CT]  (9) 

can  be  formed  to  check  if  the  system  is  completely  observable.  If  the  rank  of  Mn  equals 
n,  the  number  of  states,  then  the  system  is  completely  observable.  For  an  unobservable 
system,  the  number  of  modes  that  can’t  be  observed  is  equal  to  n-rank(M„).  For 
completeness  two  more  definitions  need  to  be  given.  A  system  is  detectable  if  all  the 
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unstable  (right  half  plane)  poles  are  observable,  and  a  system  is  stabilizable  if  all  the 
unstable  (right  half  plane)  poles  are  controllable. 

Linear  Quadratic  Regulator  f  1 6:  sect  3.3;  17:  sect  6.1] 

If  a  linear  time  invariant  system  is  completely  controllable,  the  poles  can  be  placed 
anywhere  in  the  complex  plane  using  full  state  feedback.  It  makes  sense  to  place  all  the 
poles  in  the  left  half  plane  for  stability,  but  it  is  also  desirable  to  place  the  poles  far  to 
the  left  to  make  the  transients  die  out  quickly.  By  doing  this  the  system  would  converge 
quickly  to  its  steady  state,  but  the  control  inputs  may  be  excessively  large.  The  control 
input  amplitudes  are  physically  limited  which  puts  a  constraint  on  the  speed  of  the  system 
or  how  far  left  the  poles  can  be  shifted.  The  relative  importance  of  system  speed  to 
control  input  amplitude  naturally  leads  to  an  optimization  problem. 

To  formulate  the  linear  quadratic  regulator  problem,  first  consider  the  state  space 
representation  of  a  completely  controllable  linear  time-invariant  system 

x  =  Ax  +  Bu  (1) 

with  the  full  state  feedback  law  u  =  -Kx.  Then  the  closed  loop  system 

x  =  (A  -  BK)x  (10) 

has  as  its  poles  the  eigenvalues  of 

det(XI-A  +  BK)  (11) 

which,  with  the  appropriate  feedback  gains,  can  be  shifted  to  any  desired  location.  The 
problem  is  to  place  the  poles  such  that  the  system  comes  to  its  steady  state  as  fast  as 
possible  with  the  least  control  inputs.  With  this  in  mind,  a  criteria  to  measure  the  speed 
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in  which  an  initial  state  is  returned  to  the  zero  state  (which  is  the  desired  steady  state)  is 
given  by 

|  x‘Qx  d\  (12) 

where  the  state  weighting  matrix  Q  is  a  positive-definite  symmetric  real  matrix.  The 
matrix  Q  is  chosen  to  indicate  the  relative  importance  of  each  state  component.  The 
integral  in  equation  (12)  measures  the  deviation  of  x  from  its  steady  state  value  during 
the  time  (t,-t„).  By  minimizing  this  criterion,  the  states  will  now  be  driven  to  the  steady 
state  as  fast  as  possible,  but  the  system  will  experience  large  control  amplitudes.  To 
prevent  the  large  control  amplitudes,  a  term  must  be  added  to  equation  (12)  to  include  the 
input.  Thus,  the  criterion  to  be  optimized  is  defined  by 

J  =  I'  (x'Qx  +  u'Ru)  ch  (13) 

where  the  controls  weighting  matrix  R  is  a  positive-definite  symmetric  real  matrix  which 

determines  the  relative  importance  of  the  control  inputs.  This  integral  is  called  the  linear 

quadratic  regulator  cost  function  or  performance  index.  The  problem  now  lies  in  finding 

the  input  u  that  minimizes  the  performance  index  J. 

From  optimization  theory,  the  feedback  law  u  =  -Kx  that  minimizes  equation  (13) 

as  t,  — 00  is 

K  =  R  'BTP  (14) 

The  matrix  P  is  the  positive-definite  solution  (which  is  unique)  of  the  algebraic  Riccati 
equation 

AP  +  PA  +  Q  PBR  B  P  =  0  (15) 

By  requiring  the  matrix  Q  to  be  positive  definite,  all  the  states  will  appear  in  the  xTQx 
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term.  This  ensures  that  all  initially  unstable  poles  are  observable  (which  is  a  requirement 
for  stability).  By  defining  HTH  =  Q  and  requiring  [A,H]  to  be  detectable,  the  closed  loop 
system  will  always  be  stable. 

As  was  mentioned  before,  the  designer  has  the  freedom  to  chose  the  weighting 
matrices  Q  and  R.  Equation  (15)  shows  that  the  choice  of  Q  and  R  will  influence  P,  the 
solution  of  the  algebraic  Riccati  equation.  Notice  from  equation  (14)  that  the  matrices 
P  and  R  play  a  role  in  determining  the  feedback  matrix  K.  Since  the  feedback  matrix  K 
shifts  the  poles  of  the  system,  the  elements  of  Q  and  R  become  the  design  parameters  for 
pole  placement. 

LQR  Stability  Margins  [17:  chapter  7] 

An  important  property  of  linear  quadratic  regulators  is  that  they  guarantee  robust 
stability  margins.  In  this  section,  both  single-input  single-output  (SISO)  and  multi-input 
multi-output  (MIMO)  systems  are  shown  to  have  guaranteed  margins.  But  first,  a 
relationship  that  is  simply  stated  here  but  thoroughly  derived  in  reference  [17]  is 
introduced.  By  manipulating  the  regulator  equations  (1),  (14)  and  (15),  the  relationship 
[I  +  Rl/2K(-jwI  -  A)  'BR  ,/:]T[I  +  Rl/:K(jwI  -  Ay'BRl/:]  >1  (16) 

known  as  the  Kalman  Inequality  can  be  derived. 

SISO  Margins.  For  a  SISO  system,  the  Kalman  Inequality  is 

[1  +  rl/:k(-jwl  -  A)  'br  l/2]T[  1  +  r,/:k(jwl  -  A)  'br  l/:]  >1  (17) 

where  the  lowercase  bold  letters  are  vectors.  This  can  be  simplified  to 

1 1  +  k(jwl  -  A)  'b  |  >1  (18) 

A  block  diagram  for  the  SISO  system  is  shown  in  figure  2.  The  closed  loop  transfer 
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function  of  a  SISO  system  is 

TF  =  G°"' _  (19) 

where  G(jw)=(jwl-A)'b  and  H(jw)=k.  The  denominator  of  the  transfer  function  for  this 
system  is  known  as  the  return  difference  function  and  is  given  by  l+k(jwI-A)  ‘b. 


Figure  2  LQR  SISO  Block  Diagram 


The  Kalman  Inequality  in  equation  (18)  implies  that  when  the  open  loop  transfer  function 
k(jwl  -  A) 1  b  (or  GH)  is  plotted  on  a  polar  plot,  the  magnitude  is  always  outside  a  unit 
circle  centered  at  the  (-1,0)  point.  In  other  words,  the  allowable  region  for  LQR  is 
outside  a  unit  disk  centered  at  the  (-1,0)  point.  On  the  polar  plot,  the  intersection  of  the 
open  loop  plot  and  a  unit  circle  centered  at  the  origin  gives  the  phase  margin  (PM).  The 
maximum  gain  margin  (GM)  is  defined  as  the  amount  the  gain  can  be  increased  before 
the  system  goes  unstable  (or  on  a  polar  plot,  when  it  crosses  the  (-1,0)  point).  In  LQR 
design,  the  system  will  not  go  unstable  because  the  plot  cannot  enter  the  unit  circle.  This 
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is  easier  to  see  with  a  picture.  In  figure  3  the  minimum  GM  is  at  point  A  and  the 
minimum  PM  is  measured  from  the  real  axis  to  the  lines  going  through  points  B.  The 
GM  must  always  be  greater  than  1/2,  which  is  inverse  of  the  point  where  the  polar  plot 
crosses  the  real  axis  (or  the  amount  the  gain  can  be  changed  before  the  plot  moves  from 
point  A  to  the  (-1,0)  point).  Remembering  that  LQR  is  restricted  from  entering  the  dotted 
unit  circle  shown  in  figure  3,  the  smallest  PM  possible  is  the  60“  angle  measured  to  point 
B.  This  shows  that  the  guaranteed  stability  margins  for  LQR  are 

1/2  <  GM  <  °°  (20) 

-60  '<  PM  <60“  (21) 

MIMO  Margins.  In  the  MIMO  case,  the  margins  are  first  showm  for  R=I,  and 
then  extended  to  include  any  choice  of  R.  For  a  MIMO  system  with  R=I,  the  Kalman 
Inequality  reduces  to 

[I  +  K(jwl  -  A)  'B]‘[I  +  K(jwl  -  A)-'B]  >  I  (22) 

This  relation  holds  true  if  and  only  if  the  minimum  singular  value  of  the  return  difference 

matrix  is  greater  than  or  equal  to  one.  That  is 

a  [I  +  K(jwl  -  A)  'B]  >  1  (23) 

The  gain  and  phase  margins  defined  for  the  SISO  case  do  not  extend  their  usefulness  to 

MIMO  systems.  More  meaningful  stability  concepts  called  independent  gain  and  phase 

margins  are  defined  by  Ridgely,  Banda,  and  Yeh.  They  give  the  following  definition: 

Independent  gain  margins  are  limits  within  which  the  gains  of  all  feedback 
loops  may  vary  independently  at  the  same  time  without  destabilizing  the 
system,  while  the  phase  angles  remain  at  their  nominal  values. 
Independent  phase  margins  are  limits  w  ithin  which  the  phase  angles  of  all 
loops  may  vary  independently  at  the  same  time  without  destabilizing  the 
system,  while  the  gains  remain  at  their  nominal  values  [17:3-73]. 
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Figure  3  Polar  Plot  of  LQR  Margins 
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The  independent  gain  and  phase  margins  for  a  MIMO  system  are  given  by 


_L  <  IGM  <  __ L 

1+a  1 -a 

and 


(24) 


-2-sin1  (4)  <  IPM  <  2-sin  1  (l)  (25) 

where  a  =  minimum  singular  value  of  the  return  difference.  From  equation  (23),  the 
smallest  value  of  a  is  one.  Using  a=l  in  equations  (24)  and  (25)  give,  as  in  the  S1SO 
case,  guaranteed  independent  margins  of 

1/2  <  IGM  <  «  (26) 

-60  <  IPM  <  60°  (27) 

Safonov  and  Athans  [19]  proved  that  these  margins  are  the  same  for  a  general  diagonal 
R  matrix  for  the  special  case  when  Q  >  0  and  all  perturbations  in  the  various  feedback 
loops  are  noninteracting.  For  this  case,  any  choice  of  R  could  be  used  since  a 
nondiagonal  matrix  can  always  be  replaced  with  a  diagonal  equivalent.  If  either  of  the 
two  constraints  are  violated  (which  they  typically  are),  the  only  form  for  R  that  gives  the 
guaranteed  stability  margins  is  R  =  pi  (where  p  is  any  scalar). 

The  guaranteed  stability  margin  property  of  LQR  is  very  desirable  for  the  pole 
placement  problem.  However,  the  relationship  of  the  weighting  matrices  of  LQR  to  the 
location  of  the  closed  loop  pole  is  not  well  understood.  One  method  of  placing  poles 
with  LQR  is  to  vary  Q  and  R  by  trial  and  error  and  watch  the  movement  of  the  poles, 
but  this  is  not  very  efficient.  In  the  next  section,  an  algorithm  is  presented  that  will  find 
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the  Q  and  R  matrices  that  places  the  closed  loop  poles  as  close  as  possible  to  a  set  of 
desired  poles  thus  simultaneously  obtaining  pole  placement  and  stability  robustness. 
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III.  Development  of  a  Robust  Pole  Placement  Algorithm 


The  important  properties  of  ensuring  a  system  is  stable  and  has  guaranteed 
margins  are  desirable  advantages  of  LQR.  In  this  section,  an  algorithm  is  presented 
which  finds  appropriate  LQR  weighting  matrices  to  allow  the  closed  loop  poles,  placed 
with  the  LQR  technique,  to  approach  a  set  of  desired  poles.  If  a  desired  pole  is  outside 
the  allowable  LQR  region,  the  algorithm  finds  the  pole  inside  the  LQR  region  (the 
achievable  pole)  that  is  closest  to  the  desired  pole.  The  algorithm  also  allows  a  weighting 
to  be  placed  on  the  achievable  poles  which  forces  the  higher  weighted  achievable  poles 
closer  to  the  desired  poles.  This  can  be  an  advantage  in  third  or  higher  order  systems, 
since  in  higher  order  systems  the  LQR  region  depends  on  the  location  of  each  pole.  The 
algorithm  tries  to  minimize  the  distance  between  desired  and  achievable  poles  but  might 
not  be  able  to  place  all  poles  exactly.  If  a  certain  pole  has  a  heavier  weighting,  the 
algorithm  tries  to  place  it  closer  to  the  desired  pole.  When  the  weighted  pole  is  shifted 
closer  to  its  desired  location  the  LQR  achievable  regions  change,  causing  a  shift  in  the 
other  non-weighted  poles  away  from  their  desired  locations.  Although  the  algorithm  may 
not  give  the  exact  overall  desired  pole  spectrum,  it  gives  the  designer  the  ability  to 
enforce  a  required  pole  placement.  The  algorithm  also  delivers  a  system  that  is  robust, 
stable,  and  close  to  the  required  pole  specifications. 

Defining  the  Problem 

Using  LQR  design,  the  closed  loop  system  is 
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X  =  (A  -  BK)x 


(10) 


with  the  optimal  feedback  matrix  defined  as 

K  =  R'BtP  (14) 

Using  this  feedback  the  closed  loop  system  is 

x  =  (A  -  BR  BTP)x  =  A*x  (28) 

and  the  closed  loop  (achievable)  poles  are  the  eigenvalues  of  A*.  Thus,  the  closed  loop 
eigenvalues  depend  only  on  Q  and  R  (since  P  is  a  function  of  Q  and  R). 

The  goal  of  the  algorithm  is  to  minimize  the  cost  function 

n 

/-  E  o) 

i  =  l 

where  the  pole  weights  V,  are  chosen  by  the  designer.  In  the  LQR  problem,  the  state 

weighting  matrix  is  required  to  be  positive  semidefinite  (Q  >  0)  and  the  control  weighting 

matrix  positive  definite  (R  >  0).  Q  can  be  made  positive  semidefinite  by  defining 

Q=H'H.  Likewise,  R  can  be  made  positive  definite  by  defining  R=MTM  and  requiring 

M  to  be  square  or  tall  (i.e.  more  rows)  and  have  full  column  rank.  Since  only  a  special 

case  upholds  the  guaranteed  stability  margins  for  a  diagonal  R,  let  R  be  defined  as  pi. 

As  mentioned  before,  the  system  is  always  stable  if  the  pair  [A,H]  is  detectable. 

If  Q  and  R  are  varied  by  giving  small  perturbations  to  H  and  M,  the  response  of 

d  J  *  3  J  * 

the  pole  difference  cost  function  f  can  be  determined  by  - and - A  gradient 

d  H  d  M 

^  j*  ^  j  • 

search  technique  will  find  the  minimum  value  of  f  when  -  and  -  is  smaller 

M  d  H  d  M 

than  a  preset  tolerance. 
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Algorithm 

The  Algorithm  to  place  robust  stable  poles  is  given  by  the  following  steps: 

1 )  Define  the  problem  by  specifying  A,  B,  and  the  desired  pole  locations. 

2)  Choose  the  pole  weighting  parameters  V,. 

3)  t  rovide  an  initial  guess  for  H  and  p  (since  R=pl). 

4)  Calculate  Q  and  R. 

5)  Solve  the  Riccati  equation  and  find  the  gain  matrix  K. 

6)  Calculate  the  achievable  poles. 

7)  Calculate  the  pole  difference  cost  function  f . 

8)  Starting  the  second  time  through,  check  to  see  if  the  magnitude 

of  J  anc|  \f  f  I  are  |ess  then  a  prescribed  tolerance 

0H  0M  '  K" 

(where  k  is  the  iteration  number).  If  so  stop. 

9)  Perturb  H  and  M  and  go  to  step  4. 

Using  the  Algorithm  With  MATLAB 

All  the  examples  were  run  using  MATLAB  on  a  Compaq  286  computer  or  a  VAX 
1 1/780.  MATLAB  is  a  powerful  controls  software  package  that  allows  the  user  to  write 
custom  macros  (functions)  or  "m"  files  that  can  be  integrated  with  other  MATLAB 
routines.  Six  new  "m"  files  were  written  for  the  execution  of  the  algorithm.  The  function 
FMINS  is  also  used  for  the  gradient  search.  The  six  "m"  files  are  described  below: 

1.  LQRPP  -  Linear  Quadratic  Regulator  Pole  Placement  is  the  main  body  for  the 
algorithm.  The  designer  calls  this  function  and  sends  the  A  and  B  matrices,  the  desired 
pole  placements,  and  the  pole  weighting.  It  returns  the  achievable  poles,  the  gain  matrix 
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K,  and  the  weighting  matrices  Q  and  R.  It  assigns  the  initial  values  to  H  and  p  (for  Q 
and  R),  sets  up  the  plotting  feature,  and  initializes  the  gradient  search  by  calling  FMINS. 

2.  PPFUNC  -  Pole  Placement  Function  is  the  function  that  is  minimized  by  the 
gradient  search  routine  (FMINS).  This  function  calculates  the  eigenvalue  difference  and 
includes  any  pole  weighting  specified  by  the  designer.  Care  is  taken  to  direct  the 
achievable  poles  toward  the  closest  desired  poles.  Each  achievable  pole  calculated  during 
the  gradient  search  is  plotted  by  this  function. 

3.  PPMAKEH  -  Pole  Placement  MAKE  H  is  a  function  called  by  PPFUNC  to  make 
the  symmetric  H  matrix  from  the  vector  that  is  being  perturbed. 

4.  PLOTINIT  -  Plot  Initialize  is  called  from  LQRPP  to  set  up  the  plot  axis  size  and 
plot  the  desired  poles. 

5.  PPINIT  -  Pole  Placement  Initial  is  called  by  LQRPP  to  set  up  a  vector  containing 
the  upper  triangular  values  of  an  nxn  identity  matrix.  These  values  are  sent  to  PPFUNC 
to  form  the  initial  Q  and  R  matrices. 

6.  STARTPP  -  Called  by  LQRPP  to  define  global  variables. 

Use  is  made  of  MATLAB's  gradient  search  routine  that  employs  the  Nelder-Mead 
simplex  algorithm  [19].  Also,  MATLAB's  powerful  graphing  capability  is  used  to  plot 
points  of  the  gradient  search  as  it  converges  to  the  achievable  poles.  Both  initial  guesses 
of  H  and  M  are  taken  to  be  the  identity  matrix  I.  All  the  "m"  files  used  to  run  the 
examples  in  the  next  section  are  included  in  appendix  A. 

As  the  example  problems  got  larger  than  third  order,  the  run  time  on  the  Compaq 
286  took  one  to  several  hours.  In  particular,  the  sixth  order  F-4  example  ran  for 
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seventeen  hours  before  the  PC  ran  out  of  memory.  The  VAX,  however,  was  usually  a 
factor  of  six  faster  than  the  PC;  the  F-4  example  took  only  three  hours  to  run  and  the 
third  order  example  finished  in  ten  minutes  (clock  time).  On  the  higher  order  examples 
the  plotting  feature  was  turned  off  to  speed  up  the  run  time. 
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IV.  Examples 


Five  examples  are  now  presented  to  show  the  effectiveness  of  the  robust  pole 
placement  algorithm. 

Example  l  -  First  Order  SISO  Case 

The  First  order  SISO  is  included  as  an  example  because  the  LQR  regions  and  the 
performance  index  can  be  found  in  closed  form.  The  closed  form  solutions  allow  a 
convenient  check  for  the  algorithm's  solution. 

Problem  Setup.  The  First  order  SISO  case  can  be  formulated  by  letting  the  A 
matrix  and  b  vector  be  scalars.  Then  the  open  loop  state  space  description  is 

x  =  ax  +  bu  (29) 

Using  the  full  state  feedback  u=-kx,  the  closed  loop  system  is 

x  =  (a  -  bk)x  (30) 

The  closed  loop  poles  or  eigenvalues  are 

detOvI  -  a  +  bk)  =  0  (31 ) 

or 

Jtcl  =  a-bk  (32) 

Thus,  the  solution  to  the  differential  equation  in  (30)  is  simply 

x  =  x„eA'  =  x0e“bk”  (33) 

where  x„  is  an  initial  condition  for  x.  For  the  optimal  solution,  k  is  found  from  the 
Riccati  equation;  but  k  can  also  be  found  in  a  closed  form  solution  from  the  quadratic 
performance  index.  Putting  equation  (33)  into  the  feedback  equation  gives 
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u  =  -kx0e<abkl' 

Combining  equations  (2 1 )  and  (4),  the  quadratic  performance  index  is 
J  =  x02  JH  (q  +  rk2)e2<a  bk”  dl 

=  xc2(q  +  rk2)  £  e2(“ ' bk)l  dt 

-  ^r^l[e--bki; 

2(a-bk)  1  Jo 


(34) 

(35) 

(36) 

(37) 


For  the  pole  placement  method  in  this  report,  the  control  weighting  matrix  r  is  chosen  to 
be  p.  In  this  example  some  insight  can  be  gained  by  watching  r  vary.  With  the  optimal 
k  given  by  equation  (40)  and  letting  b=l,  the  optimal  closed  loop  pole  is 


Note  that  the  closed  loop  pole  is  always  negative,  so  the  system  will  always  be  stable. 
No\*  let  q  or  r  in  equation  (41)  approach  infinity 
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O  —  oo  \  —  oo 

”  Cl  opt 

r— °°  X,  ,  —  -|al 

cl  opt 


(42) 


As  q-*°°  the  system  gets  faster  and  the  closed  loop  poles  move  further  to  the  left.  This 
demonstrates  that  q  is  indeed  the  weighting  for  the  speed  of  response  term  -  increase  q 
and  the  system  gets  faster.  On  the  other  hand,  the  value  of  r  determines  the  control 
effort.  If  control  is  expensive,  the  designer  would  make  r  big.  From  equation  (42),  as 
r— 00  the  closed  loop  pole  is  put  at  the  negative  absolute  value  of  a.  Thus,  if  the  open 
loop  system  is  stable  and  r— °°,  the  pole  does  not  move  and  no  control  is  used.  If  the 
op>en  loop  system  is  unstable,  the  optimal  solution  is  to  shift  the  pole  to  its  stable  mirror 
image. 

Equations  (40)  and  (41)  show  that  this  example  is  really  just  a  one  parameter 
system  that  depends  on  1  .  By  letting  1  =  a  and  varying  a  from  zero  to 
approximately  infinity  (a  *  00  since  r  *  0,  r=0  would  signify  unlimited  control  power), 
the  achievable  regions  for  the  closed  loop  poles  can  be  found.  By  looking  at  the  open 
loop  transfer  function 

ghol  =  k  (si  —  a)  1  b  =  hs.  (43) 

s-a 

the  open  loop  pole  is  located  at  a.  A  root  locus  for  the  LQR  achievable  poles  is  drawn 
in  figure  4  for  the  stable  and  unstable  case.  Observe  that  both  closed  loop  systems  are 
stable  (as  they  had  to  be)  and  the  locus  can  not  enter  the  area  between  zero  and  point  -a. 
That  region  is  restricted  for  all  closed  loop  poles  regardless  of  the  choice  for  Q  and  R  . 
Since  the  cost  function  that  was  minimized  was  the  LQR  performance  index,  the  only 
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allowable  closed  loop  poles  are  to  the  left  of  point  -a.  To  see  why,  look  at  the  gain  and 
phase  margins  for  the  system.  From  equation  (43),  the  margins  are  calculated  as 


f  \ 


GM  =  20  log 


opt 

a 

V  / 


PM  =  180 -tan 


v'Vp,  -a; 


where  GM  is  in  dB  and  PM  is  in  degrees.  If  a  pole  were  placed  in  the  restricted  region, 
it  would  violate  the  guaranteed  phase  margin  property  of  LQR. 

Now  that  the  allowable  regions  of  LQR  are  known,  consider  the  two  cases  of  a 
stable  and  an  unstable  open  loop  pole.  For  each  case,  let  the  desired  pole  be 


v  = a  - k 


(44) 


so  that 
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(45) 


k  =  a  -  a 


dc* 


Stable  Case.  Letting  a=-5  and  Xdcs=-7  gives 
2 


Shol  = 


s+5 


-125^-5 


=  a  -  X.  =  -5+7  =2 

dcs 


(46) 


solving  for  q  gives 

q  =  (kopi  +  5);  -  25  =  49  -25  =  24  (47) 

Then  the  achievable  pole  is 

K*  =  K,  =  -  -  724  -  25  -  -7  (48) 

The  PM  and  GM  are  guaranteed  to  be  at  least  60“  and  -6  dB.  From  the  Bode  plot  in 
Figure  5,  the  gain  margin  for  k  =  2  is  infinite.  This  result  implies  that  the  system  is  stable 
and  robust  for  any  positive  gain. 

Unstable  Case.  Letting  a=5  and  Xdes--7  gives 


gh 


ol 


12 
s  -  5 


k  =  a  -  X,  =12 

dcs 


(49) 


and 

k„P,  =  V25  *q  -  5 


(50) 


Solving  for  q  gives 

q  =  (k,r,-5)‘-  25  =  49  -25  =  24 


(51) 
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Figure  5  Bode  Plot  for  Stable  Case 

then 

K*  =  Ki.  =  =  -  V25  -  24  =  -7  (52) 

The  Bode  plots  for  this  system  are  shown  in  figure  6.  The  GM  is  now  [-7.6,  °°]  dB  and 
the  PM  is  65.4  degrees,  which  both  meet  the  guaranteed  margins  as  expected.  An 
interesting  result  in  the  unstable  case  is  when  the  system  uses  the  least  control  effort. 
That  is  when  the  unstable  pole  is  shifted  to  its  mirror  image.  Figure  7  shows  the  Bode 
plot  for  this  case,  where  Xdc,  is  now  at  -5.  The  stability  margins  are  at  [-6,  °°]  db  and  60 
degrees-  exactly  the  minimums  guaranteed  by  LQR. 

Algorithm  Results.  Using  the  results  of  the  last  section,  the  pole  difference  cost 
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Figure  6  Bode  Plot  for  Unstable  Case  k=  1 2 
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Figure  7  Bode  Plot  for  Unstable  Case  k=2a 


function  can  be  expressed  as 


J'  -  f  (53) 

where  the  pole  weighting  parameter  V  is  not  necessary  since  there  is  only  one  pole. 
Since  X.(ks  and  a  are  fixed,  f  is  a  function  of  q  only.  Taking  the  derivative  with  respect 
to  q  and  setting  it  equal  to  zero  gives  the  minimum.  Solving  for  q  gives 

q  =  *-L~a2  (54) 

If  the  desired  pole  is  in  the  achievable  region,  the  function  is  minimized  w.  r,  /  equals 
zero.  When  the  desired  pole  is  outside  the  achievable  region,  the  function  is  minimized 
when  q  equals  zero;  this  will  give  the  achievable  pole  location  at  -/af 

Both  cases  (stable  and  unstable)  were  run  on  MATLAB  using  the  "m"  files  listed 
in  appendix  A.  Figures  8  and  9  show  the  algorithm  plot  for  the  gradient  search  for  the 
best  achievable  pole.  The  process  starts  with  the  initial  value  of  q=l. 

For  all  cases  the  algorithm  successfully  places  the  pole  at  the  exact  closed  form 
location.  If  the  desired  pole  w'ere  in  the  restricted  region  between  zero  and  minus  five, 
the  achievable  pole  should  be  placed  at  negative  five.  Figure  10  show's  this  to  be  the 
case.  Table  I  compares  the  results  from  the  algorithm  with  the  closed  form  solutions. 
The  closed  form  solutions  for  the  second  unstable  case  were  calculated  the  same  as  the 
other  two  cases. 


30 


First  Order  SISO  System  Poles 

6 

4 

L.  ' 

Q  -  *  -  - - - - - - -  .  -  - - 

-4 

-6 

-8  -  -  -  •  - -  -  --  - . 

-8-6-4-2  0  2  4  6  8 

Figure  8  Algorithm  Plot  for  Stable  Case 
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Figure  10  Algorithm  Plot  for  Unstable  Case  with  Desired  Pole  at  -4 


Table  I  First  Order  SISO  Results 


Stable 

Open  Loop 
Pole 

Method 

a 

^ach 

n 

q 

Algorithm 

-5 

-7 

-7 

2 

24 

Calculated 

-5 

-7 

-7 

2 

24 

Unstable 

Open  Loop 
Pole 

Algorithm 

5 

-7 

-7 

12 

24 

Calculated 

5 

-7 

-7 

12 

24 

Unstable 

Open  Loop 
Pole 

Algorithm 

5 

-4 

-5 

10 

0 
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5 

-4 

-5 

10 

0 
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Example  2  -  Second  Order  SISO  Case 

As  the  order  of  the  problem  increases,  the  allowable  LQR  regions  become  more 
complex,  but  with  the  second  order  case  the  regions  can  still  be  calculated  in  closed  form. 
This  allows  a  check  for  the  algorithm  to  see  if  the  achievable  pole  is  indeed  the  closest 
point  in  the  LQR  region  to  the  desired  pole. 

Problem  Setup.  The  second  order  SISO  problem  can  be  formulated  by  applying 
a  force  to  a  frictionless  mass  as  depicted  in  figure  1 1.  The  equation  of  motion  for  this 
system  is 

u  =  mz  (55) 

By  letting  x,  =  z  and  x,  =  z  the  state  space  representation  is 

/- 

x. 


0  1 
0  0 


x. 


(56) 
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With  m=l,  the  system  is  defined  as 


0  1 
0  0 


"o 

1 


u  =  Ax  +  Bu 


(57) 


A  block  diagram  of  the  system  is  seen  in  figure  12,  where  the  full  state  feedback  is 


u  = 


-Kx 


(58) 


Since  u  =  z  ,  equation  (58)  can  also  be  expressed  as 


z  +  k,z  +  k,  z  =  0  (59) 

which  is  the  characteristic  equation  for  the  system.  To  see  this,  look  at  the  transfer 
function  method  for  finding  the  characteristic  equation.  For  the  system  in  figure  12,  the 
characteristic  equation  is 

1  -gh  =  0  (60) 


or 


1  *KfsI  -  A]-'B  -  0 
Expanding  equation  (61)  gives 

si  -  A 


5  0 

0 

l" 

5  -1 

0  5 

0 

0 

0  V 

si  -  A]'1 


s  1 
0  s 


fsl  -  A]  'B  = 


"  1  1  " 

s  S* 

0 

s: 

0  1 

s 

1 

s 

(61) 


(62) 


(63) 


(64) 
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kl  =  <v 


k2  =  2t>, 


(69) 


The  open  loop  system  (gh)  is  K[sl  -  A]  'B,  or 

k^s  +  k 

gh  =  - - - 

s* 

Using  k,  and  k:  from  equation  (69)  and  shifting  to  the  frequency  domain  gives 


gh(jco)  = 


2(>J(o  +  co; 


-or 


(70) 


(71) 


Dividing  through  by  wn:  to  normalize  the  problem  and  defining  w  =  —  results  in 

0) 

n 

gh  =  J  :2_\^  (72) 

-0)  “ 

Normalizing  the  problem  this  way  will  show  that  the  optimization  problem  is  not  a 
function  of  frequency  (u>n).  The  open  loop  system  can  be  represented  in  phaser  form  as 


gh(jw) 


/.(tan'1  2£w  -  180deg) 


(73) 


0) 

To  find  the  PM  for  this  system,  set  the  magnitude  equal  to  one  (|gh|=l)  and  solve  for 
the  frequency.  Choosing  the  real  solution  for  the  frequency  gives 


IT 


(74) 


Putting  this  value  into  the  angle  term  of  equation  (73)  and  subtracting  from  -180  degrees 
gives  the  PM 


PM  =  tan  '2C[2C: +  (4C4  +  irf  (75) 

By  normalizing  the  frequency,  the  system  has  been  reduced  to  a  one  parameter  problem. 
Finding  the  LQR  Regions.  On  a  Nyquist  plot,  LQR  design  prohibits  the  locus  of 
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points  of  GH(jto)  from  entering  a  unit  circle  centered  at  (-1,0).  In  figure  3,  if  the  plot  of 
GH(jto)  were  exactly  on  the  dotted  unit  circle  (centered  at  (-1,0))  LQR  would  have  its 
minimum  stability  margins.  From  vector  addition,  the  radius  of  the  dotted  unit  circle  in 
figure  3  is  1 +GH.  Thus,  a  mathematical  constraint  on  a  Nyquist  plot  from  LQR  design 
is 


Id  +  GH)|  >  1 
In  this  second  order  example. 


1  ~gh  = 


s2  +  k,s  +  k,  s2  +  2  Cwns  + 


S'  s* 

And  in  the  normalized  frequency  domain 

1  +  2  Cwj  -  (i) ' 


(1  *  gh) 


—0) 


Putting  equation  (78)  into  the  constraint  equation  (76)  leads  to 

w4  *(4C:  -  2)<jj‘  +  1  >  o) 4 
or 

(4C:  -2)w:  *  1  >0 

This  constraint  will  hold  true  for  any  frequency  only  if 
4C:  -  2  >  0 

Thus,  the  damping  ratio  must  meet 


(76) 


(77) 


(78) 


(79) 


(80) 


(81) 


C>J_=  0.7071  (82) 

The  relationship  of  £  and  con  to  poles  in  the  complex  plane  is  depicted  in  figure  13. 
Using  the  relation  cos  X  =  0,  where  0  is  measured  as  in  figure  13,  the  achievable  region 
for  LQR  design  is  where  0  <  45'.  This  result  is  plotted  in  figure  14.  Now  that  the  LQR 
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Figure  14  LQR  Barriers 
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achievable  region  is  known,  the  desired  poles  need  to  be  chosen.  Consider  the  following 
system  specifications: 

•  peak  overshoot  <  0.456 

•  rise  time  h  <  0.44  seconds 

•  settling  time  ts  <  4.6  seconds 

•  robust  system 

Then,  to  just  meet  the  system  performance  requirements,  «n  =  4.12  and  £  =  0.2425 
f20:sect  2.7.4].  From  the  diagram  in  figure  13,  the  desired  pole  is  found  to  be  at  \dts  = 
-1  +  4j.  The  desired  pole  meets  the  first  three  specfications,  but  has  a  PM  of  only  27°. 
To  find  the  best  achievable  pole  (one  that  is  robust  and  closest  to  the  desired  pole),  the 
LQR  barrier  must  be  examined.  Using  the  top  half  of  the  plot  in  figure  14,  the  equation 
of  the  LQR  barrier  is 


The  point  in  the  LQR  region  that  is  closest  to  the  desired  pole  is  at  the  intersection  of  the 
barrier  and  a  line  perpendicular  to  the  barrier  passing  through  Xdes.  Thus,  the  equation  of 
the  line  perpendicular  to  the  barrier  and  passing  through  -l+4j  is 

y  =  x  +  5  (84) 

The  intersection  of  equations  (83)  and  (84)  yields  the  achievable  pole,  or 

Kt,  =  -2-5  +  2.5j  (85) 

The  achievable  pole  gives  the  following  system  characteristics: 

•  Mr  =  0.04 
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•  4  =  0.5 1  seconds 


•  ts  =  1.84  seconds 

•  C  =  0.7071 

•  (on  =  3.54 

•  PM  =  65.5" 

The  Phase  Margin  was  calculated  with  equation  (75)  and  can  be  verified  with  the  Bode 
plot  in  figure  15.  The  Nyquist  plot  in  figure  16  shows  that  the  system  does  stay  outside 
the  unit  circle  centered  at  (-1,0). 

Notice  that  the  PM  is  not  exactly  60”.  The  LQR  design  guarantees  only  a 
minimum  PM  of  60”,  but  larger  margins  are  often  achieved.  That  is  precisely  the  case 
in  this  example.  To  show  why  LQR  can  not  produce  a  PM  of  60”  for  this  system. 


Figure  15  Bode  Plot  for  Second  Order  SISO  System 


equation  (75)  was  solved  for  £  (with  PM=60°).  The  value  £=0.612  was  found  which 
corresponds  to  the  barrier  angle  0=52.2°.  Using  the  same  analysis  of  dropping  a 
perpendicular  line  from  the  desired  pole  to  the  barrier  yields  X  =  -2.3  +  2.98j  for  a 
PM=60°.  This  gives  con=3.76  which  allows  the  gains  and  a  Nyquist  plot  to  be  calculated. 
The  Nyquist  plots  for  the  LQR  system  and  the  PM=60°  system  are  shown  in  figure  17. 
Notice  that  the  PM =60°  plot  violates  the  unit  circle  and  thus  can  not  be  generated  by 
LQR  design. 

Algorithm  Results.  The  algorithm  was  run  on  MATLAB  and  generated  the  plot 
in  figure  18.  In  the  plot,  x  marks  the  desired  pole  while  o  shows  the  achievable  pole. 
Notice  that  the  points  calculated  during  the  gradient  search  do  not  move  outside  the  45° 
lines.  Table  II  gives  the  results  for  this  example.  The  gains  for  the  closed  form  solution 
were  calculated  with  equation  (69).  The  results  show'  that  the  algorithm  is  able  to 
accurately  find  >vath.  The  MATLAB  "m"-file  used  to  solve  the  second  order  SISO  system 
is  included  in  appendix  A.  Run  time  on  the  Compaq  286  was  4  minutes. 


Table  II  Second  Order  SiSO  Results 
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Figure  16  Nyquist  Plot  for  LQR  Design  -  Second  Order  System 
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second  order  system  poles 


Figure  18  Algorithm  Plot  With  LQR  Barriers 
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Example  3  -  Third  Order  SISO  Case 

In  the  second  order  example  the  LQR  barriers  could  be  found  in  closed  form  since 
the  solution  depended  only  upon  £.  The  third  order  system  cannot  be  reduced  to  a  single 
parameter  and  solved  in  closed  form,  but  the  results  of  the  previous  sections  show  that 
the  technique  used  is  accurate.  Three  cases  are  presented  in  this  third  order  example;  the 
first  chooses  desired  poles  inside  the  LQR  region,  the  second  shows  the  effects  of 
weighting  an  actuator  pole  and  the  effects  on  the  other  closed  loop  poles,  and  the  third 
case  examines  ihe  effects  on  a  pole  that  is  both  outside  the  allowable  LQR  region  and 
close  to  the  imaginary  axis. 

Problem  Setup.  If  an  actuator  is  added  to  the  second  order  problem  it  becomes 
a  third  order  system.  The  system  can  be  modeled  as  in  figure  19.  Letting  m=l  and 
defining  the  states  as 
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Figure  19  Third  Order  System  Block  Diagram 
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X 


r 

- 

z 

z 

where  u  =  (-u  +  u.)a  and  z  = 

z 

z 

u 

J 

1  *  = 

il 

m 

the  state  space  representation  is 


0 

1  0 

0 

x  = 

0 

0  1 

x  + 

0 

0 

1 

CO 

1 

O 

a 

u.  =  fk,  k,  kj  x 


(86) 

(87) 


The  actuator  pole  is  set  by  physical  constraints,  but  for  the  sake  of  this  example  let  a=l. 


Then,  from  1  +K[sI  A]  B,  the  characteristic  equation  is 
s‘  +  ( 1  +  k,)s:  +  k2s  +  k,  =  0 


(88) 


Factoring  the  characteristic  equation  gives  three  roots,  one  from  the  actuator  and  two  from 
the  plant. 

Case  I.  For  this  case  the  desired  poles  should  be  within  the  LQR  achievable 
region.  In  the  second  order  example,  the  barriers  made  45'  lines  with  the  real  axis.  Since 
the  third  order  LQR  allowable  regions  are  unknown,  choose  a  desired  location  within  the 
second  order  system  region.  Let 


^=[-0-5+0.5j  -0.5-0.5j  -1]  (89) 

Using  the  function  PLACE  on  MATLAB,  the  closed  loop  poles  can  be  placed  in  the 
desired  locations  and  the  gain  matrix  K  calculated.  The  resulting  gains  are 


K  =  [0.5  1.5  1.0] 


(90) 


With  these  gains,  the  characteristic  equation  is 
s'  +  2s:  +  1.5s  +  0.5  =  0 


(91) 
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whose  roots  are  the  same  as  A.dcs.  From  the  Bode  plot  in  Figure  20,  the  PM=65.5”  and  the 
GM=C,°.  The  system  is  within  the  limits  guaranteed  by  LQR  design,  so  the  algorithm 
should  place  the  poles  exactly.  The  algorithm  generated  the  plot  in  figure  21  and  did 
achieve  the  desired  pole  location.  To  achieve  the  desired  poles,  the  Q  and  R  are 


R  =  8.293  and  Q  = 


2.073  0.829  -0.003 

0.829  2.067  -0.006 

-0.003  -0.006  0.0 


Case  2.  In  this  case,  the  desired  poles  are  chosen  outside  of  the  LQR  allowable 
region.  The  first  run  is  done  with  no  weighting  on  the  actuator  pole,  and  the  second  run 
weights  the  actuator  pole  three  times  as  heavily  as  the  other  poles.  The  desired  poles 
need  to  lie  outside  the  45”  second  order  system  barrier.  Thus,  let 

K,  =  f-3+5j  -3-5j  -10]  (92) 

Using  the  function  PLACE  on  MATLAB  gives 

K  =  [34  9.4  0.6]  (93) 

For  these  gains,  the  Bode  plot  in  figure  22  shows  the  PM=53.34”.  The  algorithm  was  run 
on  MATLAB  to  find  the  LQR  achievable  poles  and  plotted  in  figure  23.  The  x's  denote 
the  desired  poles,  the  dots  are  the  achievable  pole  choices  generated  during  the 
algorithm's  gradient  search,  and  the  +’s  are  the  actuator  pole  choices  from  the  algorithm. 
It  can  be  seen  that  the  desired  poles  were  outside  the  LQR  allowable  region.  The 
achievable  poles  and  phase  margin  are 

Xach  =  [-3.48  -  4.52i  -3.48  +  4.52i  -10.78]  (94) 

PM  =  62.45”  (95) 
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Figure  20  Bode  Plot  for  Desired  Poles  -  Case  1 
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Figure  21  Algorithm  Plot  -  Case  I 
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To  see  the  effects  of  weighting  a  pole,  the  same  desired  poles  were  used  but  the 
actuator  pole  was  weighted  three  times  as  heavily  as  the  plant  poles.  The  algorithm  found 
the  best  LQR  achievable  poles  as 

Xach  =  [-3.62-4.30i  -3.62 -4.30i  -10.53] 

The  gradient  search  from  the  algorithm  is  plotted  in  figure  24.  Again,  the  x’s  are  the 
desired  poles  and  the  dot  and  plus  symbols  denote  the  gradient  search  attempt  to  find  the 
best  achievable  pole.  All  achievable  pole  choices  on  the  plot  are  within  the  LQR 
allowable  region,  and  the  pattern  of  the  achievable  pole  choices  gives  an  outline  for  the 
LQR  barriers.  However,  the  LQR  regions  depend  on  the  location  of  all  the  closed  loop 
poles,  so  the  algorithm  plot  can  not  give  a  true  estimate  of  the  LQR  barriers.  Using  the 
gains  found  by  the  algorithm,  the  Bode  plot  shown  in  figure  25  was  generated. 

Table  Ill  lists  the  results  of  the  three  different  runs  for  case  two.  Notice  from 
Table  III  that  the  Phase  Margin  increased  14.5%  for  the  non-weighted  actuator  condition 
and  16.3%  for  the  weighted  actuator  condition.  The  weighted  actuator  pole  moved  2.5% 
closer  to  the  desired  location  compared  to  the  non-weighted  actuator  pole,  but  at  the 
expense  of  the  other  poles.  When  the  actuator  pole  was  weighted,  the  other  two  plant 
poles  were  1 .5%  less  in  magnitude  (con)  and  4.8%  less  in  damping  angle  (0)  from  the  non- 
weighted  locations.  There  is  a  trade  off  in  pole  location  accuracy,  and  the  designer  must 
decide  if  the  benefits  gained  from  the  weighted  pole  outweigh  the  penalties  to  the 
remaining  poles. 
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Figure  24  Algorithm  Plot  -  Case  2  (  Actuator  Weighted) 
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Figure  25  Bode  Plot  for  Achievable  Poles  -  Case  2  (Actuator  Weighted) 


Table  III  Third  Order  System  Case  2  Results 


Pole  Location 

Gains 

PM 

Desired  Pole  Location 

-3  ?  5i 
-10 

[34  9.4  0.6] 

53.34 

No  Weighting  on  Poles 

-3.48  ¥  4.52i 
-10.78 

[35.35  10.82  0.78] 

62.45 

Actuator  weighting 

-3.62  ¥  4.30i 
-10.53 

[33.29  10.79  0.78] 

63.75 

Case  3.  Like  case  2,  the  plant  poles  are  chosen  outside  the  LQR  achievable 
regions.  They  are  also  chosen  with  poor  PM  to  demonstrate  a  large  increase  in  PM  with 
a  small  shift  in  position.  The  actuator  pole  is  chosen  at  a  location  of  -2.5.  The  desired 
plant  poles  are  >wdes=[-.2+.75i  -.2-.75i].  Using  the  MATLAB  function  PLACE,  the  gain 
matrix  for  the  desired  poles  was  found.  This  allowed  a  Bode  plot  and  the  PM  to  be 
found.  The  Bode  plot  for  this  system  is  seen  in  figure  26.  Running  the  algorithm  to 
locate  the  achievable  poles  produced  the  plot  in  figure  27.  In  this  case,  the  Q  and  R 
matrices  from  the  achievable  poles  in  case  two  were  used  to  reduce  the  run  time.  The 
Bode  plot  for  the  achievable  poles  is  in  figure  28.  Using  the  function  MARGIN  on 
MATLAB  the  PM  was  found  to  be  60.74°.  This  can 

also  be  seen  from  the  Bode  plot  in  figure  28.  The  results  are  compared  in  table  IV. 
Notice  that  for  a  6.2%  penalty  in  magnitude  (con)  and  a  24.5%  smaller  damping  angle  (0) 
there  is  a  110%  increase  in  PM.  Figure  29  shows  the  location  of  the  plant  desired  and 
achievable  poles. 
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Figure  26  Bode  Plot  for  Desired  Poles  -  Case  3 
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Figure  27  Algorithm  Plot  -  Case  3 

Figure  28  Bode  Plot  for  Achievable  Poles  -  Case  3 


Table  IV  Third  Order  System  Case  3  Results 


-ii-U - — - 

Pole  Location 

PM 

-0.2  -  0.75i 
-0.2  +  0.75i 
-2.5 

28.90“ 

^ach 

-0.4  -  0.61  i 
-0.4  +  0.6  li 
-2.77 

60.74“ 
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Figure  29  Comparison  of  Desired  and  Achievable  Poles 


Example  4  -  Comparison  with  Haney  and  Stein  method  f  1 5] 

Harvey  and  Stein  present  a  procedure  to  asymptotically  place  eigenvalues  and 
eigenvectors  using  optimal  feedback.  Their  technique  is  to  formulate  the  state  weighting 
matrix  Q  such  that  the  transmission  zeros  of  the  system  are  located  at  the  desired  pole 
locations.  Then,  as  the  control  weights  tend  to  zero  (or  the  feedback  gains  increase  to 
infinity),  the  closed  loop  poles  migrate  towards  the  desired  pole  locations. 

This  technique  has  one  major  fault;  the  number  of  transmission  zeros  available  to 
be  placed  is  n-m  (#  states  -  #  control  inputs).  As  some  of  the  closed  loop  poles  migrate 
towards  the  desired  locations,  others  head  off  to  infinity.  The  poles  that  head  off  to 
infinity  happen  to  be  the  actuator  poles.  Harvey  and  Stein's  solution  to  this  problem  is 
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to  monitor  the  actuator  poles  as  the  gains  are  increased  and  choose  the  gains  that  put  the 
actuator  poles  at  their  required  locations.  They  use  an  F-4  aircraft  inner  loop  lateral  axis 
design  example  to  demonstrate  their  procedure.  The  states  and  controls  are 


x 


V, 


P 

4> 

5 

r 

5 

a 


stability  axis  roll  rate 
stability  axis  yaw  rate 
angle  of  sideslip 
bank  angle 
rudder  deflection 
aileron  deflection 


(96) 


u 


&rc  rudder  command 
6  aileron  command 


(97) 


w'hich  gives 


-.746 

.387 

-12.9 

0 

.952 

6.05 

0 

0 

.024 

-.174 

4.31 

0 

-1.76 

-.416 

0 

0 

.006 

-.999 

-.0578 

.0369 

.0092 

-.0012 

x  + 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-20 

0 

20 

0 

0 

0 

0 

0 

0 

-10  _ 

0 

10 

or 


x  =  A  x  +  B  u 


From  Mil  Spec  8785B,  the  desired  poles  are 


•  Roll  subsidence  mode  -4.0 


•  Dutch  roll  mode  -0.63  t  ?.42j 

•  Spiral  mode  -0.05 


and  the  actuators  lie  at 
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•  Rudder  actuator 


-20 


•  Aileron  actuator  -10 

Two  points  must  be  made  about  the  Harvey  and  Stein  example: 

1.  The  control  weighting  is  defined  as  pR  and  the  procedure  is  to  let  p  tend 
towards  zero.  But  to  allow  the  actuator  poles  to  achieve  the  desired  locations, 
the  procedure  has  to  be  iterated  for  many  different  values  of  p.  In  their 
example,  Harvey  and  Stein  found  that  a  value  of  p=0.0025  put  the  actuator 
poles  close  to  the  desired  locations.  As  p  goes  to  zero  the  plant  poles  go  to 
the  desired  locations.  But  the  actuator  pole  position  determines  the  value  of 
p  and  hence  the  closed  loop  plant  pole  locations. 

2.  Since  the  actuator  poles  head  out  to  infinity  as  the  gains  are  increased  (or  p 
is  decreased),  Haney  and  Stein  modified  the  initial  pole  locations  in  the  A 
matrix.  The  open  loop  actuator  poles  were  moved  to  (-10,  -5)  to  allow 
freedom  to  increase  the  gains  before  the  actuator  poles  moved  past  the  desired 
locations.  There  has  to  be  some  play  in  the  gains  to  allow'  the  other  poles  to 
migrate  towards  the  transmission  zeros,  and  Harvey  and  Stein  modified  their 
example  to  start  the  actuator  poles  to  the  right  of  their  desired  locations. 

When  the  F-4  data  w'as  run  with  the  algorithm  presented  in  this  report,  the  A 
matrix  was  changed  back  to  the  original  actuator  values  (-20,  -10)  to  reflect  the  correct 
locations  for  the  actuator  poles.  This  example  was  run  on  the  VAX  using  MATLAB. 
Both  the  algorithm  and  Harvey  and  Stein  were  able  to  place  all  the  poles  close  to  their 
desired  locations.  Figures  30  and  3!  show  the  desired  and  achievable  pole  locations  for 
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both  the  algorithm  results  and  Harvey  and  Stein’s  solution.  Table  V  compares  the  results 
for  this  example.  Notice  that  Harvey  and  Stein  achieved  poles  close  to  the  desired  plant 
locations,  but  the  algorithm  was  12%  closer  to  the  desired  aileron  actuator  pole  and  4% 
closer  to  the  desired  rudder  actuator  pole.  The  robustness  for  a  MIMO  system  can  be 
measured  using  singular  values  as  discussed  in  chapter  II.  Both  Harvey  and  Stein’s 
solution  and  the  algorithm’s  solution  were  examined  for  robustness  using  singular  values, 
and  the  singular  value  plots  for  these  two  cases  are  shown  in  figure  32.  Two  robustness 
tests  can  be  checked  using  singular  values.  By  defining  K  fsl-A]  'B=T'  (where  the  * 
depends  on  which  case  is  being  tested),  the  tests  are 


a[I  *T*]  >0  dB 
a[I  +  T"1]  >  -6  dB 


(98) 


These  first  test  is  plotted  in  figures  33  and  34,  and  the  second  test  is  plotted  in  figures 
35  and  36. 


Table  V  Example  Four  Results 


Harvey  and  Stein 

Algorithm 

-3.810 

-3.998 

-0.049 

-0.091 

-0.727  +  2.358j 

-0.669  *  2.365j 

-10.43 

-10.025 

-22.44 

-20.053 

m 

f.  132  .882  -1.576  -.026  -.681  .026] 

[-.306  -1.389  .729  .039  .107  -.089] 

■ 

[-.524  -.42  2.827  -.021  .013  -.86  j 

[  .409  .858  -.060  .035  -.044  .239J 

PM 

[-53”,  53°] 

[-60°,  601 

CiM 

[-5.56  dB,  19.7  dB] 

[-6.02  dB.  oo  dB] 

57 


40 

30 

singular  value  plot  -  Algorithm(-)  Harvey  and  Steini ..) 

20 

10 

- 

£ 

0 

- 

E 

5j 

-10-  - 

-  ■  - 

- 

-20 

- 

-30 

- 

-40 

50 

-60 

10  ■’ 

10-  10  1  10°  10'  10- 

103 

w  -  rad/sec 

Figure  32  Singular  Value  Plots  for  Example  4 


Using  the  minimum  singular  values  from  figures  33  and  34,  the  independent  gain  and 
phase  margins  were  calculated  with  equations  (24)  and  (25). 

As  was  pointed  out  in  chapter  II,  R  can  be  diagonal  for  guaranteed  stability 
margins  only  for  a  special  case.  That  happens  when  the  disturbance  matrix  is  diagonal 
or  each  one  is  independent  of  the  others.  If  the  disturbances  are  not  all  independent,  the 
restriction  is  even  tighter;  R  must  equal  pi  for  guaranteed  margins.  In  this  example, 
p=25.38.  When  R  is  allowed  to  be  diagonal,  where  the  diagonal  elements  do  not  equal 
each  other,  the  guaranteed  robustness  breaks  down.  This  proves  that  the  disturbances  in 
the  various  feedback  loops  interact  with  each  other  [19].  A  plot  of  robustness  test  1  for 
the  diagonal  R  matrix 
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Figure  34  Minimum  Singular  Value  Test  1  for  Algorithm 
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Figure  36  Minimum  Singular  Value  Test  2  for  Algorithm 
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Figure  37  Singular  Value  Test  for  a  Diagonal  R 


Figure  38  Singular  Value  Test  for  a  Symmetric  R 


R  = 


0.2424 

0 


0 

0.8630 


is  given  in  figure  37.  The  IPM  is  reduced  to  (-50.4°,  50.4°),  well  below  the  guaranteed 
minimum  of  (-60°,  60°).  However,  there  is  a  benefit  for  the  reduction  in  stability  margins 
-  the  achievable  poles  were  placed  exactly  in  the  desired  locations.  Using  a  diagonal  R 
also  gives  one  more  parameter  to  design  with.  Let’s  go  one  step  further.  What  if  R  was 
allowed  to  be  symmetric?  This  would  give  three  variables  in  R  that  could  be  changed, 
thus  adding  two  extra  degrees  of  freedom  to  the  original  problem.  The  algorithm  was 
modified  to  let  R  be  symmetric  and  returned 


R 


0.0805  0.0581 
0.0581  0.0563 


The  first  robustness  test  in  equation  (98)  is  plotted  in  figure  38  for  this  symmetric  R 
matrix.  All  the  poles  were  placed  in  the  desired  location,  but  the  margins  were  reduced 
even  further,  the  IPM  =  (-40”,  40”)  and  the  !GM  =  (-4.5  dB,  9.97  dB).  If  the  designer 
had  some  freedom  to  relax  the  stability  margins,  a  choice  of  a  diagonal  or  symmetric  R 
matrix  could  be  used  to  place  all  the  poles  in  their  desired  locations.  But  some  care  must 
be  taken  if  R*pl  If  any  one  value  in  R  is  much  lower  than  the  others,  the  corresponding 
parameter  in  the  control  vector  would  have  unacceptably  large  input  amplitudes  [7:124], 
The  cross  coupling  also  has  constraints.  If  Z  represents  unidirectional  cross  coupling, 
then  for  a  diagonal  R 
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This  is  mentioned  to  highlight  some  of  the  potential  problems  that  occur  when  R*pl.  In 
this  example  all  the  elements  in  R  are  of  the  same  order,  but  each  design  case  is  different. 
The  pole  placement  accuracy  could  outweigh  the  problems  encountered  with  a  diagonal 
or  symmetric  R. 

The  Harvey  and  Stein  technique  placed  the  zeros  in  the  desired  pole  locations  by 
fixing  Q  and  varied  R  to  move  the  poles.  The  algorithm  allowed  both  Q  and  R  to  vary, 
and  achieved  better  results  in  both  pole  placement  and  robustness.  Harvey  and  Stein’s 
results  shown  in  table  V  were  calculated  for  their  modified  A  matrix.  If  their  gains  were 
applied  to  the  actual  aircraft  A  matrix,  their  actuator  pole  locations  would  be  further  left 
and  the  stability  margins  would  be  worse.  This  example  again  demonstrates  the  pole 
placement  accuracy  using  the  method  presented  in  this  report.  In  the  next  example,  an 
aircraft  w;th  open  loop  level  II  flying  qualities  is  examined  to  see  if  the  algorithm  can 
improve  flying  qualities  with  a  robust  controller. 
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Example  5  -  Aircraft  Flying  Qualities  Improvement 

This  example  demonstrates  the  algorithm’s  ability  to  improve  the  flying  qualities 
of  an  aircraft.  The  aircraft  in  this  example  is  a  fictitious  model  similar  to  an  A-4D,  with 


characteristics  primarily  obtained  from  reference  [21],  Appendix  B  lists  the  stability 
derivatives  used  in  this  example.  Considering  the  longitudinal  axis  with  only  the  elevator 


deflection  as  input,  the  state  space  representation  is 


-0.0129 

-3.7292 

0 

-32.2 

1l 

-0.0002 

-0.8167 

0.9984 

0 

0 

x  = 

-0.0003 

-1.6903 

0.0563 

0 

X  + 

1.56 

0 

0 

1 

0 

0 

where 


u 


velocity 


x 


a 

q 

e 


angle  of  attack 
pitch  rate 

pitch  angle 


(99) 


and 

u=5t.  (elevator  deflection) 

The  open  loop  poles  for  this  system  lie  at 

-3.84  +  1.1 9j  (short  period) 

-0.0026  *  0.0397  lj  (phugoid) 


(100) 


which  give 


Cr  =  0.307 
Cph  -  0.065 


0) 


n  sp 


1.25 


“n  Ph 


=  0.04 
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This  flight  profile  (cruising  flight)  is  classified  as  category  B  by  the  Mil-Spec  users  guide 
[22: 14].  For  category  B,  a  level  I  flying  qualities  rating  is  defined  when  0.30<£sp<2.0  and 
Cph>0.04.  The  short  period  frequency  requirements  are  given  in  figure  39.  The  regions 
of  flying  qualities  levels  shown  in  figure  39  are  for  n/a,  where  n  is  the  load  factor 
[23:538-540].  Roskam  [23]  gives  an  approximation  to  n/a  as 
n  ~Z 

_  -  _ 2  (102) 

«  g 

where  Za  =  Zw.u.  With  the  data  in  appendix  B,  n/a  is  found  to  be  16. 1 1  g’s/rad.  Plotting 
o)n  vs  16. 1 1  in  figure  39  falls  within  the  level  II  region  (point  A).  The  values  for  £sp  and 
Cph  fail  within  the  level  I  range. 

The  goal  is  to  make  the  system  both  level  I  and  robust  with  full  state  feedback. 
Thus,  the  desired  poles  were  chosen  to  give  level  I  damping  and  place  the  frequency 
requirement  in  the  center  of  the  level  I  region.  Point  B  in  figure  39  corresponds  to  the 
desired  short  period  frequency.  The  desired  poles  and  associated  £  and  «n  are 
- 1 . 1 2  ?  3.50j  CsP  =  0-305  a)D  sp  =  3.67 

-0.0056  +  -0.073j  CPh  =  0.076  o>nph  =  0.073 

Using  PLACE  on  MATLAB,  the  unique  gains  were  calculated  that  gave  the  Bode  plot 
for  the  desired  poles  in  figure  40.  The  gain  and  phase  margins  from  figure  45  are 
GM  =  °o  and  PM  =  38°.  The  GM  is  fine,  but  the  PM  is  not  quite  high  enough.  The 
algorithm  can  not  achieve  the  desired  poles  but  it  should  get  close.  When  this  system 
was  run  on  MATLAB  with  the  algorithm  m-files,  the  achievable  poles  were  found  to  be 
-2.096  *  2.389j  (short  period) 

-0.215  *  0.043j  (phugoid) 
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Figure  39  Short  Period  Frequency  Requirements  -  Category  B  [23:541] 


6: 


The  natural  frequency  and  damping  are 


CsP  =0.6595  tonsp  =  3.179 
Cph  =  0.98  conph  =  0.219 

which  are  all  within  the  level  1  region  and  close  to  the  desired  values.  In  figure  39,  point 
C  shows  the  short  period  frequency  location  for  the  achievable  poles.  The  achievable 
system  Bode  plot  shown  in  figure  41  gives  the  GM=°°  and  the  PM =68".  Thus,  using  the 
algorithm  for  pole  placement  gave  a  robust  level  I  system. 
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Figure  40  Bode  Plot  for  Desired  Poles  for  Example  5  Aircraft 
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Figure  41  Bode  Plot  for  Achievable  Poles  for  Example  5  Aircraft 


V.  Conclusions  and  Recommendations 


This  report  has  examined  the  pole  placement  properties  of  linear  quadratic 
regulators  as  determined  by  the  weighting  matrices  of  the  performance  index.  A  n2w 
eigenvalue  difference  cost-function  to  minimize  the  difference  between  sets  of  achievable 
poles  and  desired  poles  was  introduced.  This  led  to  the  formulation  of  an  algorithm  that 
employed  a  gradient  search  to  optimize  the  eigenvalue  difference  performance-criterion. 
This  algorithm  gives  a  powerful  method  for  selecting  weighting  matrices  for  pole 
placement  designs. 

Five  examples  were  used  to  demonstrate  the  algorithm.  The  use  of  simple  systems 
in  the  First  two  examples  allowed  closed  form  solutions  which  compared  exactly  to  the 
solutions  of  the  algorithm.  The  third  order  example  showed  the  effects  of  weighting  a 
desired  pole.  The  heavier  weighted  pole  was  forced  closer  to  the  desired  pole  at  the 
expense  of  the  less  weighted  poles.  A  method  presented  by  Harvey  and  Stein  was 
compared  to  the  algorithm  given  in  this  report,  and  the  results  show  that  the  algorithm 
accurately  provided  robust  eigenvalues  for  a  MIMO  system.  In  the  last  example,  the 
algorithm  was  used  to  improve  aircraft  flying  qualities.  These  examples  showed  that  the 
algorithm  can  simultaneously  place  desired  poles  and  achieve  stability  robustness. 

This  algorithm  only  employed  an  eigenvalue  optimization  performance  criterion. 
The  effects  of  optimizing  both  an  eigenvalue  and  an  eigenvector  difference  on  achievable 
pole  location  should  be  investigated.  The  full  effects  of  the  LQR  weighting  matrices  on 
pole  locations  is  still  not  well  understood.  This  algorithm  could  be  modified  to  perturb 
selected  values  in  the  weighting  matrices  to  see  the  effects  on  the  closed  loop  poles.  One 
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improvement  that  could  be  made  is  in  the  gradient  search  routine.  The  algorithm  simply 
uses  the  Nelder-Mead  method  on  MATLAB  which  calculates  gradients  numerically,  but 
a  more  efficient  method  could  be  found.  A  possible  way  to  increase  the  efficiency  is  to 
use  the  eigenvalue  difference  weighting  parameter  to  automatically  weight  the  poles  with 
the  largest  eigenvalue  difference.  The  m-files  could  undoubtedly  be  programmed  more 
efficiently,  and  a  Fortran  source  code  would  be  useful  to  users  who  do  not  have  access 
to  MATLAB.  This  method  is  applicable  to  many  flight  control  system  design  problems 
and  could  be  used  to  improve  the  robustness  properties  of  existing  systems. 
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Appendix  A.  MATLAB  m-files 


Appendix  A  includes  listings  of  the  "m"  files  used  on  MATLAB  to  run  the  examples 
presented  in  this  report.  The  following  MATLAB  functions  were  also  used  but  are  not 
included  in  this  appendix: 

FMINS  -  The  gradient  search  minimization  function  in  the  MATLAB  files. 

PLACE  -  Used  to  find  the  gains  required  to  place  a  pole  at  a  desired  location.  This 
function  is  in  the  Control  System  Toolbox  files. 

MARGIN  -  This  function  finds  the  gain  and  phase  margins  for  a  given  system.  It  is  in 
the  Control  System  Toolbox  files. 

LOGSPACE  -  A  function  in  the  MATLAB  files  used  to  generate  a  logarithmically 
spaced  frequency  vector. 

BODE  -  Located  in  the  Control  System  Toolbox  file  and  used  to  generate  Bode  plots. 


function  [ea,k,Q,R|=LQRPP(A,B,ed,V) 

%  LQRPP  Pole  placement  using  Linear  Quadratic  Regulator  design. 

7c  fea,k]=LQRPP(a.b,ed)  will  attempt  to  place  poles  in  a 

7c  desiied  location  (ed).  If  the  desired  location  is  outside 

7c  the  allowable  LQR  region,  LQRPP  minimizes  the  distance 

7c  between  the  desired  poles  (ed)  and  the  achievable  poles  (ea). 

7c  A  desired  pole  can  be  weighted  heavier  to  indicate  its 

7c  relative  importance.  LQRPP  will  give  a  weighted  pole 

%  priority  in  the  pole  placement  routine.  The  technique  is  to  use 

7c  a  gradient  search  to  minimize  the  difference  between  the 

7c  desired  and  LQR  achievable  poles  for  the  system 


72 


x  =  Ax+Bu  with  full  state  feedback  u  =  -Kx 


% 

% 

%  The  required  inputs  are 
%  A  -  the  state  space  A  matrix 

%  B  -  the  state  space  B  matrix 

%  ed  -  a  vector  containing  the  desired  pole  locations 
% 

7c  There  is  one  optional  input 
%  V  -  the  pole  weighting  matrix  (diagonal) 

% 

%  There  are  four  output  arguments  (the  last  two  are  optional) 

% 

7c  ea  -  a  vector  containing  the  achievable  poles 
%  K  -  the  feedback  gain  matrix  to  get  the  achievable  poles 

%  Q  -  the  LQR  state  weighting  matrix 

7C  R  -  the  LQR  control  weighting  matrix 

% 

startpp 

a_s=A; 

b_s=B; 

ed_s=ed; 

[n,m]=size(A); 

(l,p]=size(B); 

%  assign  initial  values  to  v,  q,  and  r 
if  nargin  >  3,  v_s=V; 

else,  v_s=cye(n);  end 
h=ppinit(n); 

m=l;  7c  initial  value  for  rho 

x=[h  m]: 

plotinit(n.ed_s) 

H=fmins('ppfunc\x,le-2); 
for  i=  1  :n 

plot(real(ea_s(i)),imag(ea_s(i)),’og’) 

end 

title( 'Gradient  Search  For  Achievable  Poles’) 

xlabelfreal') 

ylabel(’imag’) 

pause 

k=k_s; 

ea=ea_s; 

Q=q_s; 

R=r_s; 
hold  off 
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function  J=ppfunc(x) 


% 

% 

% 


used  by  LQRPP  to  Find  the  value  of  the  eigenvalue  difference  cost  function 
This  function  calculates  f  and  is  called  by  the  gradient  search  (FMINS)  routine. 


[s,t]=size(x); 

fn,ma]=size(a_s); 

fnb,mb]=size(b_s); 

(nd,md]=size(ed_s); 

xh=x(  1  :(s-l),l); 

xm=x(s,l); 

m=xm*eye(mb); 

h=ppmakeh(n,xh); 

q_s=h’*h; 

r_s=m'*m; 

k_s=lqr2(a_s,b_s,q_s,r_s); 

ea_s=sort(eig(a_s-b_s*k_s)); 

%find  the  poles  that  are  closest  to  each  other  and  take  their  difference 

eatmp=ea_s; 

edtmp=sort(ed_s); 

i=l; 

J=0; 

while  i  <  n 
fnn,mm]=size(eatmp); 
z=abs(eatmp(  1  ),lrones(nn,mm)-edtmp); 
fvt.is]=min(z); 

7c  Find  if  there  are  any  weighted  poles 
err=fed_s-edtmp(is)*ones(nd,md)); 
fvp,ip]=min(abs(err)); 

J=J+vtA2*v_s(ip,ip); 
eatmp=eatmp(2:nn); 
if  is==l,  edtmp=edtmp(2:nn); 

else  if  is==nn,  edtmp=edtmp(l:(nn-l)); 
else,  edtmp=[edtmp(  1  :is- 1  );edtmp(is+ 1  :nn)]; 
end;end 
i-i+1; 
end 

err=(ed_s-edtmp(  1  )*ones(nd,md)); 
f  vp.ip]  =min(abs(err)); 

J=J+(abs(eatmp-edtmp))A2*v_s(ip,ip); 
for  i=  1  :n 

plot(real(ea_s(i)),imag(ea_s(i)),\r') 

end 
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startpp 


/C 
% 

%  Used  by  LQRPP  to  initialize  the  global  variables  used  in  the  function 


j=sqrt(-l); 

global  a_s  b_s  k_s  q_s  r_s  ea_s  ed_s  v_s 


function  h=ppmakeh(n,x) 

% 

%  PPMAKEH  is  used  by  LQRPP  to  make  the  symmetric  h  matrix  from  a  vector 
%  n  is  the  size  of  the  matrix  and  x  is  a  vector  that  contains  the  upper  triangular 
%  elements  of  the  h  matrix 
s=l; 

for  i=l:n; 
for  j=i:n; 
h(i,j)=x(s); 
h(j,i)=x(s); 
s=s+ 1; 
end 
end 


function  plotinit(n,ed) 

%  used  with  LQRPP  to  set  up  the  plotting  feature 

axis(  ’square’) 
e=sort(ed); 

if  abs(real(e(n)))  >  abs(imag(e(n))),  axisize=ceil(abs(real(e(n))))+2; 

else  axisize=ceil(abs(imag(e(n))))+2; 
end 

axis(f-axisize  axisize  -axisize  axisize]) 

plot(fO  0], [-axisize  axisize], ’-b\ [-axisize  axisize], [0  0],’-b') 

hold  on 

for  i=  1  :n 

plot(real(e(i)),imag(e(i)),’xg,) 

end 
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function  [fl=ppinit(n) 


%  used  with  LQRPP  to  initialize  the  weighting  matrices 

% 

h=n; 

x=l; 

for  i=l:n 
f(x)=l; 
for  j=  l:(h-l) 
x=x+ 1; 
f(x)=0; 

end  %j  loop 
h=h-l; 
x=x+ 1 ; 
end  %i  loop 
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Appendix  B.  Data  for  the  A-4D  Aircraft 


The  following  data  is  for  an  A-4D  flying  in  level  flight  at  15,000  ft  at  mach  0.6. 


u,  (ft/sec) 

- 

634.2 

W  (lb) 

- 

17,578 

Tu  (1/sec) 

- 

0.000225 

X„  (1/sec) 

- 

-0.01288 

X„  (1/sec) 

- 

-0.00588 

X^  f(ft/sec~)/rad] 

- 

0 

Zu  ( 1/sec) 

- 

-0.1012 

Z„  j,„  (-) 

- 

-0.001616 

Z„  (1/sec) 

- 

-0.818 

Zv  [( ft/sec  :)/rad) 

- 

0 

Mu  (1/sec-ft) 

- 

-0.000407 

d/ft) 

- 

-0.000556 

M„  (l/sec-ft) 

- 

-0.003 

Mq  (1/sec) 

- 

58.1 

\1N  (l/sec:) 

- 

1.56 

The  dynamic  equations  of  motion  are  given  on  the  following  page. 
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The  equations  of  motion  for  the  longitudinal  axis  are  [23:172;  22:54] 

<i  -u,Xw  +  g0  =  Xbb 
(u,  ~  U,ZW)6  -  Zuu  -  u,Zwa  -  (u,  +Zq)q  =  Zfi5 
q  -  M  u  +  u,  M  a  +  u,M  a  +  Mq  =  M..  6 

0  =  q 
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